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Abstract 

High-frequency (HF) stimulation has been shown to block conduction in excitable cells including neurons and cardiac 
myocytes. However, the precise mechanisms underlying conduction block are unclear. Using a multi-scale method, the 
influence of HF stimulation is investigated in the simplified FitzhHugh-Nagumo and biophysically-detailed Hodgkin-Huxley 
models. In both models, HF stimulation alters the amplitude and frequency of repetitive firing in response to a constant 
applied current and increases the threshold to evoke a single action potential in response to a brief applied current pulse. 
Further, the excitable cells cannot evoke a single action potential or fire repetitively above critical values for the HF 
stimulation amplitude. Analytical expressions for the critical values and thresholds are determined in the FitzHugh-Nagumo 
model. In the Hodgkin-Huxley model, it is shown that HF stimulation alters the dynamics of ionic current gating, shifting the 
steady-state activation, inactivation, and time constant curves, suggesting several possible mechanisms for conduction 
block. Finally, we demonstrate that HF stimulation of a network of neurons reduces the electrical activity firing rate, 
increases network synchronization, and for a sufficiently large HF stimulation, leads to complete electrical quiescence. In this 
study, we demonstrate a novel approach to investigate HF stimulation in biophysically-detailed ionic models of excitable 
cells, demonstrate possible mechanisms for HF stimulation conduction block in neurons, and provide insight into the 
influence of HF stimulation on neural networks. 
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Introduction 

Electrical signaling is fundamental to the physiological function 
of excitable cells such as neurons and cardiac myocytes. Irregular 
electrical patterns in the brain and heart can lead to life- 
threatening conditions including epileptic seizures and ventricular 
fibrillation. External stimulation can terminate these irregular 
rhythms [1,2]; however large strength stimuli are often associated 
with detrimental effects such as pain [3] and impaired cardiac 
function following defibrillation [4]. 

In the 1960s, it was shown that Idlohertz-range high frequency 
(HF) sinusoidal stimulation could reversibly block conduction in 
neurons [5]. The use of 1-40 kHz HF- induced neural conduction 
block has recendy been exploited in clinical studies for diagnostic 
and therapeutic purposes, improving bladder function [6,7] and 
mitigating pain associated with peripheral nerve activity [8-10]. 
Despite the clinical usage of HF stimulation treatment, the 
mechanisms underlying therapeutic success in these physiological 
and pathological settings are unclear. Simulation studies in 
neurons have suggested two mechanisms: reduced sodium channel 
availability due to transmembrane potential depolarization and 
persistent activation of potassium channels [8-14]. However, the 
relative significance of the two mechanisms varies with the 
properties of the neuron, as well as the specific species and model. 
Further, in simulation studies, the transmembrane potential, ionic 
currents, and channel gating variables oscillate on the fast time 
scale of the HF stimulus, varying throughout the HF stimulation 



period, such that distinguishing the precise influence of the HF 
stimulus is difficult. 

Alternatively, one can apply a multi-scale method, separating 
the fast time scale dynamics — due to the HF stimulus — and the 
slow dynamics of the excitable cell, and derive an averaged model, 
which accounts for the HF stimulus but does not contain a high- 
frequency term [15]. Using this type of approach, several studies 
have analyzed the influence of a HF stimulus in the simple 
FitzHugh-Nagumo (FHN) model [16]. Cubero and colleagues 
demonstrated that the model cell cannot repetitively fire when the 
HF stimulus amplitude-frequency ratio is above a critical value 
[17]. Ratas and Pyragas showed that this ratio also influenced 
conduction speed in a nerve axon and above a critical value led to 
conduction block [18,19]. The FHN model is minimalistic, 
reproducing many important aspects of cellular excitability [20], 
and ideal for analysis, as the model only contains two variables, 
permitting the use of standard nonlinear dynamics techniques such 
as phase-plane analysis. However, in general biophysically-detailed 
models of excitable cells are more complex than represented by 
the simple two-dimensional FHN model. 

In this study, we first illustrate the multi-scale method to derive 
the averaged FHN (AFHN) model equations and use phase-plane 
analysis to determine critical HF stimulus thresholds above which 
the model cannot exhibit repetitive firing or elicit a single action 
potential. We then extend this approach to simulate the dynamics 
of the classical Hodgkin-Huxley (HH) neuron model [21] and 
illustrate similarities and differences between the AFHN and 
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averaged Hodgkin-Huxley (AHH) model. Further, we demon- 
strate that HF stimulation alters the ionic current activation and 
inactivation dynamics, illustrating possible mechanisms for con- 
duction block in a single neuron. Finally, we simulate HF 
stimulation in a network of neurons and demonstrate that HF 
stimulation alters network synchrony and, above a critical 
stimulation strength, terminates persistent network activity, 
suggesting implications for clinical therapy. 

Results 

Averaged Fitzhugh-Nagumo model 

We begin by deriving and analyzing a simplified model of 
spiking which accounts for the influence of HF stimulation. We 
consider the FitzHugh-Nagumo (FHN) model [16] with the 
addition of a constant current / and a time-varying HF stimulus 



V = fiy) — w + I + pco cos {on) 



= £(v + /J — yw) 



(la) 



(lb) 



where ft) is the frequency of the HF stimulation, p is the HF 
amplitude-frequency ratio. 



and the dot indicates differentiation with respect to time t. The 
HF stimulation term is defined in terms of p and w with the 
foresight that the influence of the HF stimulation depends on p, 
not specifically the amplitude pw. The FHN model is a simplified 
model that reproduces many important properties and dynamics 
of excitable cells. The simplicity of the model permits a geometric 
illustration — through phase plane analysis — of many important 
biophysical phenomena such as repetitive spiking and depolariza- 
tion block. In the model, v is the dimensionless transmembrane 
potential, and ii' represents the degree of refractoriness. Through- 
out the paper, we frx e = 0.008, /? = 0.8, and y = 0.5, such that in 
the absence of any external stimuli, the neuron is excitable. 

If the period of the fast HF stimulus is much smaller than all 
characteristic times of the FHN model, according to the method of 
averaging [15], an approximation to the slow system can be 
obtained by averaging over the period of the HF stimulus. As 
shown in the Methods, the variables of the averaged Fitzhugh- 
Nagumo model (AFHN), v and w, ait governed by the following 
system of equations: 



=f(v) — w-\-I = F(v,w; pj) 



w = £(v + /? — yw) = G(v,w) 



(2a) 



(2b) 



= c(p)v- 



and 



c(p) = l- 



The AFHN mode! is very similar to the FHN model, with the 
only difference being the modification of the cubic function / that 
influences the dynamics of v. In the absence of HF stimulation, i.e., 
p = 0, the two models are identical. In the following sections, we 
investigate how the HF stimulus parameter p influences the 
properties of repetitive action potential firing. In Figure 1 , we plot 
V from simulations of the FHN model (black) for various values of 
the HF stimulation frequency ft) and compare with v from a 
simulation of the AFHN neuron model (red), for fixed values for p 
and /. In general, as O) increases, v from AFHN model becomes a 
better approximation of the average value of v from the FHN 
model simulation, validating our formulation. 

Repetitive firing in the AFHN model. In the parameter 
region considered in this study, the cell is excitable, that is, in the 
absence of an external stimulus, the cell is at rest, and the addition 
of a stimulus can induce a single or multiple action potentials. In 
this study, we wiU consider two types of applied current stimuh: a 
constant apphed current and a brief applied current pulse, in 
addition to the HF stimulation. 

We first consider the case of a constant applied current /. In 
Figure 2 A, we plot v for a constant current 7=1.3 and various 
values of p. For no HF stimulus [p = 0), the cell fires repetitively. 
Increasing the amplitude of the HF stimulation parameter p 
decreases the action potential amplitude and increases the firing 
frequency. Consistent with previous studies [17], increasing p 
further results in cessation of repetitive firing, following a single 
action potential at the stimulus onset. Conditions for cessation of 
firing are derived as follows. 

For the parameters chosen, the AFHN model has a single 
steady-state (v,,W',), which satisfies the implicit expression 



F(%,{%+ti)/r,p,i)=Q, 



(3a) 



ffl = 0.1 



- CO = 10 100 



where 



/(v) = 



271 



f[v + p sin (B)]de 



Figure 1. Validation of the AFHN model. Simulated v traces from 
the FHN (black) and AFHN (red) models for varying HF stimulus 
frequency co. Parameters: Radial frequency co is identified in each panel, 
7 = 1.3, p=l. 

doi:1 0.1 371/journal.pone.0081 402.g001 
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V V p I 



Figure 2. Repetitive firing in the AFHN model. (A) Simulated v traces for / = 1.3 and different values for p. The dashed lines indicate v = 0. (B) 
Phase-plane portrait for variable / and p. In each panel, the v-nullcline (green) is shown for 3 values of /. The w-nullcline (blue) is independent of / 
and p. (C) I-p parameter space, denoting regions of rest, repetitive firing, and block. The limit cycle lower and upper limits (/+, Eq. 6 ) and rheobase 
Uri,, Eq. 9 ) as functions of p. (D) Frequency and amplitude of action potentials, as functions of / and p. 
doi:10.1371/journal.pone.0081402.g002 



w. = (v,+/?)/y, (3b) 

and shown in Figure 3. A.s p increa.se.s, the resting potential v, 
becomes more depolarized and approaches 0 for large p. The 
degree of refractoriness also increases as p increases, such that vv* 
approaches /?/)' for large p. 

Using standard techniques from linear stability analysis [22], the 
stability of the steady-state (v,,vp,) can be determined by 
linearizing around (Vt,Wt), and evaluating the matrix of partial 
derivatives, the Jacobian J, at the steady-state. 




0 2 4 6 8 10 

P 



Figure 3. Steady-state of the AFHN model. The steady-state 
transmembrane potential and degree of refractoriness vv, are shown 
as functions of the HF stimulation parameter p. Critical values of p for 
repetitive firing p, and for evoking a single action potential following a 
brief applied current pulse p,, are identified. See text for description of 
critical values. 

doi:1 0.1 371/journal.pone.0081 402.g003 



/(V„H',) = 



When the steady-state becomes unstable, specifically the real 
part of the eigenvalues of /, c(p) — — £y > 0, a stable limit cycle 
emerges, which can be interpreted biophysically as repetitive 
action potential firing. The critical parameter value at which the 
limit cycle emerges is known as a Hopf bifurcation. Many previous 
studies have shown that in the FHN model (i.e., p = 0), as the 
applied current / increases, there are two critical values for /, /_ 
and /+ , which correspond to the onset and offset of the stable Kmit 
cycle, respectively [23-25]. Below /_, the steady-state is stable 
corresponding to the cell at rest, between /_ and /+ the steady- 
state is unstable and the cell repetitively fires, and above /+, the 
steady-state is stable again and the cell is in depolarization block 
[23]. 

In Figure 2B, we plot the nuUclines of the AFHN model for 
several values of / and p. The v-nuUcline (green) — given by the set 
of all points {v,w) such that v = 0 — is a cubic function of v, while 
the vP-nuUcine (blue) — similarly defined as the set of all points (v,w) 
such that vP = 0 — is linear, and the nuUclines intersection denotes 
the location of the steady-state. For a given value of p, increasing / 
shifts the v-nuUcline upwards, while the w-nuUcline is independent 
of both / and p. 

If / is such that the steady-state is located on the middle branch 
of the v-nuUcline, and if iv is sufficiendy slow compared to v, that 
is, e« 1, then it can be shown that the steady-state is unstable, and 
a stable limit cycle exists [23]. From a geometric illustration, we 
can anticipate a critical value of p, p,, above which a stable limit 
cycle and repetitive firing cannot exist, consistent with Figure 2A 
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(bottom panel). As p increases, the slope of the middle branch of 
the v-nuUcline decreases, and the "knees" of the nuUcline move 
towards the steady-state (Vj,^,). When the slope at (v^jW,) equals 
0, the middle branch of the nuUcline no longer exists and, 
therefore regardless of /, a stable limit cycle also does not exist. 
Using the slope of the v-nuUcline alone as a criterion for the critical 
value of p, p. 

From linear stability analysis, we can more precisely determine 
the necessary condition for a limit cycle, c(p) > ey, such that p^ is 
given by 



p, = ^/2{]^). 



(5) 



For all values of p>Ps the steady-state is always stable, 
regardless of/, as previously shown by [17]. Further, for p<p„ 
the critical stimulus upper and lower limits, /+ and 
respectively, are given by 



P>Pc 



y V V, 



(7) 



the AFHN model cannot be excited by a brief applied current, 
where is defined in Eq. 27. Using the parameters used in this 
study, 1.302. For p> p^, regardless of the magnitude of the 
stimulus pulse /, v(t) relaxes back towards the steady-state value v, 
following the applied current pulse, without a large amplitude 
excursion typical of an action potential (Figure 4A, right panel). 

For p < p^, and the strength-duration curve is approximated by 



f{%)[v*-m 
" 1 - oxp\f'(%)d] ■ 



(8) 



where the prime indicates differentiation with respect to v, such 
that 



i±{p)=Pv-'±[y 



-c(p)](c(p)-eyf^±-(c{p)-eyf\ (6) 



The /+ curves separate the regions of rest, repetitive firing, and 
depolarization block in the I-p parameter space and coalesce 
when p = Ps at a double Hopf bifurcation (Figure 2C). For the 
parameters used in this study, p,*; 1.411. 

In the regime for repetitive firing, we derive an approximation 
for the action potential frequency and amplitude in the AFHN 
model (see Methods). For a given value of p, the frequency first 
increases then decreases as / increases (Figure 2D), while the 
amplitude is constant, consistent with a relaxation oscillator. 
Increasing p increases frequency and decreases the action 
potential amplitude, consistent with Figure 2A. 

Excitability in the AFHN model. We next consider the 
excitability of the AFHN model following a brief applied current, 
in the presence of HF stimulation, by determining the .strength- 
duration curve, the relationship between the duration d of an applied 
current pulse and the minimum amplitude /q such that an action 
potential fires [26]. 

With the system initial at rest, i.e., v(f = 0) = v,, we make the 
assumption that an action potential is fired when v(/) reaches some 
threshold ^. Although it has been shown that the FHN model does 
not strictly exhibit aU-or-none threshold behavior [23], when vP is 
sufficiently slow compared with v, the middle root of the v- 
nuUcline is a reasonable approximation for an action potential 
threshold, which we show increases as p increases (see Methods for 
details and references on firing threshold, Eq. 44, Figure 4C). This 
threshold-like behavior is illustrated in Figure 4. We plot v as a 
function of time following brief ii? = 0. 1 current pulses for p = 0 and 
1. For both values of p, an action potential is elicited if v(0>^ 
during the brief pulse, while if v(?)< S during the current pulse, v 
returns to rest v,. Increasing p increases both the v threshold for 
evoking an action potential, ^, and the stimulus threshold / 
necessary to elicit an action potential (Figure 4C). 

In the Methods section, we show a critical value for p, p^., exists, 
which for all values of p. 



/'(V,) = c(p)-v2 



We plot the strength-duration curves in Figure 4B for several 
values of p. For all values of p, /o decreases linearly with d when 
presented on a logarithmic scale and approaches a constant value 
for long a relationship typical of excitable cells. For a given 
stimulus duration d, the strength required to elicit an action 
potential Iq increases as p increases. Two important values are 
typically determined from the strength-duration curves: rheobase 
(/r/,), defined as /o for an infinite duration pulse, and chronaxie 
(Tf), defined as the pulse duration having a threshold that is twice 
the rheobase. From Eq. 8, /,•/, and Tc are given by 



/./,(p)=/'(v.)[v.-«p)] 



(9) 



and 



(p) = 



ln(2) 

/'(v.)' 



(10) 



respectively. We plot Iri, and Tc as a function of p in Figure 4C. 
Both Ii-i, and are fairly constant for small p. 7^/, increases and Tj. 
decreases, as p further increases towards p^.. We also plot 7^, in 
Figure 2C for comparison with I , and note that for all values of 
p, Ii-i, < I , that is, a smaller 7 is required to elicit a single action 
potential than to ehcit repetitive spiking, as expected. We note that 
the derivation of Eq. 8 assumes the stimulus 7o is brief — that is, Eq. 
8 is strictly valid for small d — therefore, p^ should not be 
interpreted as a critical p above which no action potentials can be 
elicited by longer duration stimuli. Indeed, p,, < p,, and therefore, 
the cell can repetitively fire during long duration stimuli for 
p^<p< p,, and a single action potential can be ehcited by long 
duration stimuli for p> p^. Further, since rheobase is defined as a 
stimulus threshold for infinite d, Eqs. 9 and 10 should be 
interpreted as approximations derived from Eq. 8, which 
nonetheless provide quahtative relationships between the 
strength-duration curve parameters 7^/, and and HF stimulation 
parameter p that can be compared with a biophysically-detailed 
model, as discussed in the next section. 
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Figure 4. Excitability in the AFHN model. (A) Simulated v traces during brief d = 0.1 stimuli pulses of amplitude / for p = 0, 1, and 1.5. In 
simulations that v exceeds the threshold c, an action potential is elicited. Inset shows an expanded time course. (B) Strength-duration curve ( Eq. 8 ) 
for several values of p. (C) Rheobase (Iri,, Eq. 9 ) and chronaxie (i,., Eq. 10 ) as functions of p. 
doi:1 0.1 371 /journal.pone.0081 402.g004 



In summary, increasing the HF stimulation parameter p 
increases the thresholds for both repetitive firing and a single 
action potential, /_ and Iri,, respectively. We derive expressions 
for critical values of p and determine the influence of HF 
stimulation on the resting potential, firing frequency and 
amplitude, action potential threshold, rheobase, and chronaxie. 
These theoretical relationships provide references that can be 
compared to results from a more realistic neuron model described 
in the next section. 

Averaged Hodgkin-Huxley model 

We next derive and analyze the influence of HF stimulation on 
a biophysicaUy-detaUed model of the neuron, utilizing the 
techniques described in the previous section. We consider the 
classical space-clamped Hodgkin-Huxley (HH) neuron model of 
the giant squid axon [2 1], with the addition of an applied current / 
and HF stimulus, given by the following system of equations: 

V =[/ + pco COS (cot) - gNain^Kv -Ena)- gxn^iv - Ek) 

11a) 

-?L(v-£i)l/C 



variable, and n the potassium activation variable. Current 
conductances, reversal potentials, and gating variable dynamics 
are described in the Methods. 

Assuming that the period of the fast HF stimulation is much 
shorter than the characteristic times of the dynamics of v and the 
gating variables, as in the previous section, we approximate the 
dynamics of the slow variables by averaging over the period of the 
HF stimulus. The variables of the averaged Hodgkin-Huxley 
(AHH) model, v, m, h, and n, are governed by the following system 
of equations: 

^=[1 - gNam^hiv - Enc) - gKn*(v - Ek) - giiv- El)]/ C (12a) 
m = a,„(\ —m) — li„m = (mrj-j —m)/zi„ (12b) 



h = a.,,{\-h)- Pi,h = (ha, - h)j%i, (12c) 

n = a„(l -«)-j8„« = (wx, -«)/t„ (12d) 
where 



m = a„,( 1 - m) - P„,m = (moo - (lib) 



h = oc,,{l-h)-pi,h = ihc.-h)/Ti, 



(11c) 



av(v,p) = 



271 



oi^[v + p sin (6) /C\de, 



(12e) 



« = a„(l —«) — /?„« = («,x. —n)/z„ 



(lid) 



where HF stimulus parameters p and co are defmed as before. In 
the HH neuron model, v = V,„ — represents the transmem- 
brane voltage V„, relative to the resting potential Vycsi, m the 
sodium activation gating variable, h the sodium inactivation gating 



271 



li,[v + p&m(B)/C]de, (12f) 
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(12g 



and 



(12h) 



for xe{m,h,n}. Because of the simplicity of the FHN model, we 
could derive analytical expressions for the dynamics of the AFHN 
model variables. In the HH model, the expressions for the and 
P-^ terms that govern the dynamics of the gating variables are 
complex, and as such, it is not possible to derive analytical 
expressions for Eqs. 12e and 12f without using approximations for 
the exponential function. Therefore, Eqs. 12e and 12f are 
computed by numerical integration for particular values of v and 
P- 

As with the FHN model, we plot V,,, from simulations of the 
HH model (black) for various values of the HF stimulation 
frequency and compare with V„, = v+ Vn-st from a simulation of 
the AHH neuron model (red), for fixed values for p and / (Figure 
5). Below an HF stimulus frequency co of 5 kHz, there is significant 
disagreement between the averaged and original model. As a> 
increases, Vm from the AHH model becomes a better approxi- 
mation of the average value of K„, from the HH model simulation, 
validating the use of the averaging method. 

Repetitive firing in the AHH model. As in the previous 
section, we consider the influence of an applied current / in the 
AHH model, in the presence of HF stimulation. In Figure 6A, we 
plot Vm for different values of p, such that the neuron is 
repetitively firing, i.e., I>I . For sufficientiy large p, the neuron 
does not repetitively fire. 

We plot the I-p parameter space for the AHH model in Figure 
6B. The parameter space is qualitatively similar to the AFHN 
model, such that the range of / for which the neuron repetitively 
fires becomes smaller as p increases, and above a critical value of 
p, Pj, the neuron does not repetitively fire. In the HH model, it has 
been shown that action potential frequency increases and the 
action potential amplitude decreases for increasing / [23], and we 
find that this is true for a given value of p. For a given /, as p 
increases, in agreement with the AFHN model, action potential 
amplitude decreases. However, in contrast with the AFHN model, 
the frequency decreases as p increases (Figure 6C). 



Excitability in the AHH model. We next consider 
excitability in the AHH model following brief applied current 
pulses. Here, we consider both positive (cathodal) and negative 
(anodal) applied current stimuli. As with the FHN model, the HH 
model is known to not exhibit a strict aU-or-none firing threshold. 
However, especially for brief (0. 1 ms) pulses, the HH model 
demonstrates a threshold-like response. In Figure 7 A, we plot 
for different values of p and /. Consistent with the AFHN model, 
the V„, threshold for evoking an action potential, ^, increases for 
increasing p (left, middle panels). Further, above a critical value of 
p, Pc, an action potential cannot be evoked, regardless of / (right 
panel). Although V,,, reaches levels near 0 mV, these responses 
should not be considered action potentials, as the depolarization of 
V,„ does not arise as a consequence of the regenerative activation 
of inward currents but rather solely as a perturbation due to the 
large applied stimulus. Specifically, above p^, regardless of 
stimulus amplitude /, Vm is maximally depolarized as the end of 
the stimulus pulse and does not become further depolarized 
following the pulse. 

We also consider the influence of HF stimulation on excitability 
following anodal break stimulation, also known as post-inhibitory 
rebound. In the classical HH model (i.e., p = 0), a negative (anodal) 
applied current pulse / hyperpolarizes the steady-state resting 
transmembrane potential F,„ (Figure 7B), which permits sodium 
inactivation recovery, i.e., h moves closer to 1. Following the pulse 
offset (break), Vm returns towards the more depolarized initial 
resting potential, and due to the slower sodium inactivation 
kinetics, V„, rebound can be sufficiently large to evoke an action 
potential. As with cathodal stimulation, the threshold for 
stimulation, (determined as the magnitude of the hyperpolar- 
ization necessary for a post-inhibitary rebound), increases for 
increasing p (left, middle panels), and above a critical value of p, 
Pg an action potential cannot be evoked (right panel). C is larger 
for anodal stimulation (Figure 7D top panel, red), compared with 
cathodal stimulation (black), and the difierence increases as p 
increases, meaning a relatively larger anodal stimulation is 
necessary to evoke an action potential. Consistent with this 
finding, we find that p„ « 0.21 <p^ 0.69 pA/cm^-Hz. 

Strength-duration curves for cathodal and anodal stimulation in 
the AHH model are shown in Figure 7C. Consistent with the 
AFHN model, for a given duration, the necessary cathodal applied 
current strength /q increases as p increases. Further, rheobase Irf, 
increases and chronaxie t,. decreases as p increases, as in the 
AFHN model (Figure 7D, black traces, middle and bottom panels). 
As p approaches p^, the strength-duration curve becomes flatter, 
consistent with a decreasing chronaxie, and illustrating that for 



100 Hz 



-50 
-100 



10 kHz 

20 kHz 10 ms 



Figure 5. Validation of the AHH model. Simulated V„, traces from tfie HH (black) and AHH (red) models for varying HF stimulus frequency co. 



Parameters: (a = liif rad/s (where / is the frequency identified in each panel), / 
doi:1 0.1 371/journal.pone.0081 402.g005 



30 ^lA/cm , p = 0.1 ^lA/cm Hz. 
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0 0.05 0.1 0.15 0 50 100 150 

P (^lA/cm^Hz) /(|jA/cm2) 



Figure 6. Repetitive firing in tlie AHH model. (A) Simulated V,„ traces for / = 50 ^(A/cm^ and different values for p. The dashed lines indicate 
F„, = 0 mV. (B) l-p parameter space, denoting regions of rest, repetitive firing, and block. The limit cycle lovj/er and upper limits (/+) and rheobase 
(Irh) as functions of p. (C) Action potential frequency and amplitude, as functions of / and p. p in units of ;iA/cm" Hz. 
doi:10.1371/journal.pone.0081402.g006 



large p the magnitude of the apphed pulse, and not the duration, 
determine whether an action potential is evoked. For a given p, 
anodal rheobase is slightly larger compared with cathodal 
rheobase (Figure 7C, D). As p approaches p^, in contrast with 
cathodal strength-duration curves, the anodal curves become 
steeper, such that chronaxie increases as p increases (Figure 7D, 
bottom panel), illustrating that short duration anodal pulses 
become relatively less efTective for evoking post-inhibitory rebound 
action potentials. 

Dynamics of the AHH model. For the AFHN model, we 
demonstrate that p^ and p^ can be approximated via theoretical 
analysis of the two-dimensional dynamical system, based primarily 
on analysis of the influence of p on the phase plane. Various 
approaches have been used to simplify the HH model to a FHN- 
like two-dimensional system, often assuming fast sodium activation 
and a linear relationship between gating variables h and n for a 
given / [20]. However, we found that a similar phase plane 
analysis using this type of reduction of the AHH model was only 
moderately successful at reproducing AHH dynamics, likely due to 
the complex relationship between the gating variables dynamics 
over a wide range of / and p (not shown). 

In the AFHN model, the HF stimulation parameter p influences 
the dynamics of v through the cubic function /. In contrast, in the 
AHH model the dynamics of v are altered indirectly through the 
influence of p on the gating variables. In Figure 8A, we plot the 
steady-state activation, inactivation, and time constant curves as 
functions of F„, for different values of p. As p increases, the sodium 
activation Woo and inactivation Aco steady-state curves are shifted 
to the right, the potassium activation n^j steady-state curve is 
shifted to the left, and all three curves are less steep (Figure 8A). 
The time constants T,„, T/,, and T„ all decrease as p increases. 



Shifts in the activation, inactivation, and time constant curves 
alter the AHH system steady-state (Figure 8B). As p increases, the 
steady-state transmembrane potential Km,* becomes more hyper- 
polarized, reaching a minimum of ~ 1 0 mV hyperpolarized below the 
baseline resting potential Vrest = — 80 mV. As p increases further, 
is gradual depolarized, approaching a maximum value of 
~ 1 0 mV depolarized above Vyesi ■ The steady-state sodium activation 
gate w, decreases and approaches 0 as p increases. Despite ¥,„,* 
becoming more hyperpolarized for small p, the steady-state 
sodium inactivation gate A, also decreases, and then increases and 
approaches 1 for large p. The steady-state potassium activation 
gate n» is also complex, first increasing then decreasing and 
approaching 0 as p increases. 

Mechanisms of conduction block 

The influence of the HF stimulus parameter p on the dynamics 
of the gating variables provides significant insight into the 
mechanism of conduction block in neurons and the various 
thresholds for repetitive firing and excitability (Pj, p^, and p^) 
(Figure 8). For smaU p<Pj (the critical value of p for repetitive 
firing), the neuron can repetitively fire for some range of the 
applied current /6[/_,/+]. As p increases, the resting sodium 
channel activation gate w» and inactivation gate A, decrease and 
the resting potassium channel activation gate «, increases, which 
aU drive the neuron towards being less prone to firing. 

As p increases for small p < Pa (the critical value of p for anodal 
stimulation), the time constant for sodium channel inactivation t| 
also decreases, that is, foUowing a negative apphed current pulse, h 
will return to its resting value in a shorter amount of time. 
Combined with a more hyperpolarized resting transmembrane 
potential Vm,* and decreasing h^, the threshold for anodal 
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d (ms) d (ms) P(MA/cm2Hz) 



Figure 7. Excitability in the AHH model. (A) V„, traces following brief (0.1 ms) cathodal and (B) anodal stimulus pulses, for different values of p. 
Threshold c indicated in each panel. (C) Cathodal and anodal strength-duration curves for different values of p. (D) Cathodal and anodal threshold c 
(for 0.1 ms stimuli), rheobase, and chronaxie, as functions of p. Current pulse amplitudes in (A): 64-66 (left); 633-635 (middle); 600, 800, 1000 (right); in 
(B): 198-200 (left); 397-399 (middle); 400, 600, 800 (right); in /lA/cm^. 
doi:10.1371/journal.pone.0081402.g007 



excitation ^ increases dramatically a.s p increases (Figure 7). p = Pi, 
when is at a minimum and tJJ ha.s decreased by approximately a 
factor of 2. 

As p increases for p<Pc (the critical value of p for cathodal 
stimulation), the sodium activation curve Woo is right shifted, and 
combined with a more hyperpolarized F^,* and decreasing T^, 
results in an increasing threshold for cathodal simulation t (Figure 
7). p = Pc when V„,^f « — 90 mV and « 0 are at a minimum, 
At « 1 and «0.6 are at a maximum and gating dynamics are 
fast, that is, the time constants for sodium activation T*j and 
inactivation T;, are near 0. Rapid and persistent activation of the 
potassium current opposes the sodium current and prevents 
sufficient depolarization to reach the threshold for evoking an 
action potential. As p increases for large p > p^., rit and T* 
decrease, as gradually transitions from hyperpolarized to 

depolarized, relative to Vi-^si. For very large p»p^, all three time 
constants are essentially equal to 0, such that the gating variable 
kinetics can be defined by instantaneous functions of v. In this 
regime, the four-dimensional AHH model (Eqs. 12a-12h) is 
reasonably approximated by a one-dimensional system, for which 
a large amplitude excursion typical of an action potential is no 
longer possible. Indeed, for a system with a single stable steady- 
state (as is the case for large p), all perturbations from the steady- 
state are followed by a gradual relaxation back to rest, as observed 
for large p in Figxire 7A and B (right panels). 

The AHH model suggests that there are different mechanisms 
of conduction block, depending on the strength of the HF 



stimulus. For small p<p„ repetitive firing ceases due to decreased 
sodium channel activation (decreased m) and availability (de- 
creased h) and increased potassium channel activation (increased 
n). For intermediate ps<p<pc, gating variable dynamics are fast 
(i.e., the time constants approach 0), and therefore eliciting a single 
action potential via anodal and cathodal excitation fails due to 
rapid sodium current inactivation, in addition to decreased sodium 
activation and increased potassium activation. For large p»p^, 
sodium and potassium currents are persistently de-activated (i.e., 
mxO and nxO), preventing an action potential from being 
evoked, and V,„,* is depolarized due to the inffuence of the leak 
current. 

In summary, the inffuence of the HF stimulation parameter p 
on the properties of action potential firing in the AHH model is 
similar to that demonstrated in the AFHN model, however with 
differences in the inffuence on the resting potential and firing 
frequency. Further, simulation and analysis of the biophysicaUy- 
detailed model provides insight into the mechanisms of conduction 
block. We determined three critical values for p, which above the 
neuron cannot repetitively fire (p,), and an action potential cannot 
be evoked by cathodal (p^) or anodal (p„) stimulation. Below the 
critical values, we demonstrated that the thresholds for evoking 
repetitive firing or a single action potential increase as p increases. 

HF stimulation of an AHH neuronal network 

Finally, we consider the inffuence of HF stimulation on a 
random network of 100 neurons, each with, on average, 10 
connections, coupled via both excitatory and inhibitory synapses 
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9 P 

Figure 8. Steady-state of the AHH model. (A) Steady-state gating variables fha^, Aoo, "oo and time constants t„„ t/„ t„ as functions of K,„ in the 
AHH neuron model for different p values. (B) Steady-state values for the transmembrane potential V,„ and the gating variables (left), and gating 
variable time constants at K,,,,* (right), as functions of p. Vertical dashed lines indicate p,,, p„, and p,, (see text for description). In the top panels, the 
horizontal dashed line indicates K,,,, = — 80 mV for p = 0. Time constants in units of ms, and p in units of pA/cm'-Hz. 
doi:1 0.1 371 /journal.pone.0081 402.g008 



(see Methods for description of synaptic currents and network 
architecture). Following a single initial applied current pulse, in the 
absence of HF stimulation (i.e., p = 0), a rastergram .shows that 
most neurons in the network fire repetitively (Figure 9A), while a 
few neurons remain quiescent due to the absence of incoming 
excitatory synaptic connections (arrows). The pseudo-electroen- 
cephalogram (pEEG, Eq. 32 ) becomes disorganized after an initial 
time period of synchronization following electrical actix-ity 
initiation. 

We plot the pEEG and the corresponding power spectrum for 
the same network as p increases (Figure 9B). In general, the 
average neuron firing rate tends to decrease as p increase, 
although though this trend is not strictly monotonic with p (Figure 
9B, middle panels). In general, decreased firing rate is associated 
with increased network synchronization (Figure 9B, right panels), 
illustrated by the narrowing of the dominant peak in the power 
spectrum and increase in the synchrony measure x ( Eq. 33 ). 
When p increases further above a critical value of p, p„, network 
electrical activity ceases immediately following the initial initiation 
(Figure 9B, bottom panel, p„ = 0.05/(A/cm^-Hz for this example). 

The critical value p„ for complete cessation of network electrical 
activity was highly dependent on the specific architecture of the 
network and the relative proportion of excitatory and inhibitory 
synaptic connections (see Methods). We determined p„ as 
functions of Pexcit, the probability that each synaptic connection 
was excitatory, and the specific network architecture. For each 
value of Pcxcit, Nsii„ = 10 different networks were randomly 



constructed with the same average connection properties but 
different connections. The mean value for p„, plus/minus one and 
two standard errors of the mean (SEM, standard deviation divided 
by y/N~^,) are shown as functions of Pexcit (Figure 10). For all 
network architectures, p„ = 0 for all networks with pexcit = 0 or 
Pexcif > 0.5, that is, the network does not exhibit persistent activity, 
even in the absence of HF stimulation. In this specific type of 
network, persistent network activity requires both excitatory and 
inhibitory synaptic connections, and indeed that a majority of 
synapses are inhibitory. For 0 0.5, p„ is an inverted U- 
shaped function oipexcu, with a maximum near pexcit = 0.2. In all 
networks, the values of p„ are less than the single cell critical values 
for repetitive spiking p,, and cathodal and anodal stimulation, p^. 
and p„, respectively. 

Discussion 

Summary of main findings 

In summary, we find that HF stimulation alters the dynamics of 
excitable cells and networks, resulting in conduction block and 
electrical quiescence. In a simplified excitable cell model, we 
identify analytical expressions for HF stimulation strength critical 
values for repetitive firing and evoking action potentials. In a 
biophysically-detaUed neuronal model, we demonstrate that HF 
stimulation alters the dynamics of ionic current gating, leading to 
reduced cellular excitability and conduction block. HF stimulation 
of a neural network reduces the overall network activity and 
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Figure 9. Electrical activity in a network of AHH model neurons. (A) Rastergram of action potentials and the pseudo-electroencephalogram 
(pEEG). Arrows indicated quiescent neurons. Parameters: p = 0, = 0.3. (B) pEEG (left, Eq. 32 ), the corresponding power spectrum (middle, value 
indicates the average neuron firing rate), and synchrony measure / (right, Eq. 33 ), as functions of p. p in units of fiA/cm Hz. 
doi:10.1371/journal.pone.0081402.g009 



increases network synchronization, leading to network quiescence 
for a sufficiently large HF stimulus. 

Relation to prior work 

Previous studies have investigated the mechanisms underlying 
conduction block in neurons [8-14]. In these studies, two primary 
mechanisms for conduction block are posed: persistent potassium 
current activation of potassium opposing sodium current and 
preventing the neuron transmembrane potential from reaching a 
threshold; and reduced sodium channel availability due to a 
baseline depolarization of the transmembrane potential. It is noted 




Pexcit 



Figure 10. HF stimulation of a network of AHH model neurons. 

The critical value p„ for repetitive activity in a AHH model neural 
network, as a function of the probability of excitatory synaptic 
connections The mean p„ (thick black) over W,,,,, = 10 simulations, 
+ 1 (solid black) and 2 (thin black) SEM (standard error of the mean, 
standard deviation of i/A',,„,) are shown. Single cell critical values of p 
{Pc,Paj 3nd p,) are identified (see text for description), p in units of 
/iA/cm^-Hz. 

doi:1 0.1 371/journal.pone.0081 402.g01 0 



in several studies that these mechanisms are not mutually exclusive 
and indeed both likely play a role, depending on the species and 
specific cell type. Our findings are generally consistent with these 
previous studies, in that we find persistent potassium current 
activation and reduced sodium channel availability (or de- 
inactivation). However, a multi-scale method approach permits 
us to identify how changes in ionic current gating result in different 
critical thresholds. Further, we are able identify the influence of 
HF stimulation on the gating variable kinetics, i.e. the gating 
variable steady-state activation, inactivation and time constant 
curves, as a significant and novel factor regulating conduction 
block. 

The method of averaging has been previously used to 
investigate the influence of HF stimulation in the FHN model 
[17-19]. Cubero and colleagues previously demonstrated that HF 
stimulation can lead to cessation of repetitive firing above a critical 
threshold (a finding reproduced in the present study) [17]. Ratas 
and Pyragas determined conditions for which HF stimulation 
results in slowed and failed propagation. [18,19]. Here, we extend 
the approach of these prior studies to demonstrate the influence of 
HF stimulation on the threshold for evoking a single action 
potential and additionally to investigate HF stimulation in a 
biophysically-detaUed neuronal single-cell model and network. 
The FHN model is ideal for mechanistic studies, as the two- 
dimensional model enables phase-plane analysis and often permits 
analytical expressions for critical values. We demonstrate that 
many aspects of HF stimulation of FHN model neurons are 
similarly reproduced in the HH model, specifically qualitatively 
similar I-p parameter space for repetitive firing; influence of HF 
stimulation on the action potential threshold; and the existence of 
critical HF stimulation amplitudes above which neurons cannot 
repetitively firing or trigger a single action potential. 

However, some important properties of HF stimulation of the 
AHH model are not qualitatively reproduced by the AFHN 
model, which paints a simpler picture for the mechanism of 
conduction block. In the AFHN model, the resting potential is 
gradually depolarized and the gating variable w gradually 
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increases as p increases, suggesting that conduction block occurs to 
do the collective influence of transmembrane potential depolar- 
ization and a larger degree of refractoriness. In the AHH model, 
the resting potential is hyperpolarized for small p and depolarized 
as p increases. Additionally, the gating variable steady-state values 
and time constants are altered in a complex manner, such that 
conduction block is due to both altered refractoriness and the 
time-dependent dynamics of the refractory variables. The AFHN 
model also not does reproduce the influence of HF stimulation on 
firing frequency. In the AFHN model, the frequency increases as p 
increases. However, in the AHH model, frequency decreases, 
likely due to the reduced sodium channel availability, i.e. the 
decrease in h„, that occurs as p increases for small p<Ps- 
Understanding how HF stimulation influences firing frequency is 
significant and necessary for optimizing HF stimulation therapy, as 
frequency plays a significant role in neural computing [27]. 

Physiological significance of findings 

The term, high frequency stimulation, is often used in different 
clinical settings with different meanings. HF stimulation frequen- 
cies range several orders of magnitude from 1 00 Hz — typical of 
studies of deep brain stimulation to treat movement disorders such 
as Parkinson's disease and other neurological disorders such as 
epilepsy [28,29] — to 40 kHz — including clinical applications such 
as pain mitigation and improved bladder voiding [8-10]. An 
inherent assumption in the derivation of the averaged excitable 
cell models is that the time scale of the HF stimulus is significantly 
shorter than the time scale of cellular dynamics. We show that the 
averaged and original HH model begin to agree when the HF 
stimulation frequency is near 5 kHz (Figure 5), consistent with 
~0.5 ms (~2 kHz) time-scale for sodium channel activation, and 
there is greater agreement as the HF frequency increases. This 
suggests that the method of averaging approach may not be strictly 
appropriate to the investigation of deep brain stimulation using 
lower frequency HF in the 100-200 Hz range but highly relevant 
to the study of peripheral nerve stimulation and clinical 
applications typically utilizing kUohertz-range HF stimulation. 
Recently, Weinberg and colleagues demonstrated that HF 
stimulation in the 100-200 Hz range could block electrical 
conduction in cardiac tissue, a novel approach to terminate 
arrhythmias [30,31]. Since the time scale of cardiac dynamics is 
generally slower than neuronal dynamics, future work is necessary 
to determine the validity of the averaging method for investigation 
of the influence of HF stimulation in cardiac tissue. 

In this study, we found that HF stimulation could prevent 
persistent network electrical activity at lower HF stimulation 
amplitudes necessary for single cell quiescence, which suggests 
network activity may cease due to HF stimulation sufficiently 
reducing excitability in a subset of neurons that prevents re- 
excitation throughout the network. We only considered network 
activity that persists (or ceases) following a single initial applied 
current at the simulation onset. The network response to a 
constant, repetitive, or random (Poissonian) applied current, in 
addition to the HF stimulation, may be significantly different. We 
additionally only consider sparsely connected random network 
architectures. In this study, we found that the specific network 
architecture was highly important in determining the response to 
HF stimulation, and thus it is reasonable to speculate that the 
network response in highly connected and/or directed neural 
networks could be different from our findings. Several studies of 
models including more detailed network architecture and specific 
cell types have suggested mechanisms underlying deep brain 
stimulation treatment using HF stimulation in the 100-200 Hz 
range. Rubin and Terman demonstrated that HF stimulation of 



the subthalamic nucleus can regularize globus pallidus firing and 
eliminate pathological thalamic rhythmicity [32]. Using a systems 
theoretic approach, Agarwal and Sarma demonstrate that HF 
deep brain stimulation improves reliability of thalamic relay [33]. 
As noted above, the averaging method may not be strictiy 
appropriate in this frequency range. However, investigation of 
more physiological neural network architectures and cell types 
may suggest alternative deep brain stimulation therapies within a 
higher frequency (kHohertz) regime or provide insight into the role 
of specific network components with different responses to HF 
stimulation, as in the aforementioned studies. 

In this study, we investigate a simple network with a random 
architecture and consider the influence of HF stimulation as 
function of the relative fraction of excitatory synaptic connections. 
We illustrate a general approach to study HF stimulation in a large 
neural network which does not require simulation of the HF 
stimulation term and thus does not require a prohibitively small 
simulation time step. Here, we consider HF stimulation in the 
context of only a few network parameters. However, neural 
networks can exhibit rich and complex dynamics, and much work 
has demonstrated that the local network architecture can have 
significant influence on global behavior, e.g., the small-world 
phenomenon [34—36], and network architecture will likely 
significandy influence our findings. Additional work is necessary 
to understand the influence of HF stimulation in the context of 
networks of varying degrees of connectivity and structure. 

The critical values for repetitive firing and evoking action 
potentials are defined in terms of the HF stimulus frequency-to- 
amplitude ratio; that is, as the HF stimulus frequency increases, so 
must the HF stimulus amplitude for the same response. In a 
therapeutic device, it is ideal to minimize the amplitude of an 
applied stimulus, to minimize power consumption and mitigate 
safety issues for both the patient and device. Thus, determining the 
optimal frequency regime to minimize the amplitude for optimal 
HF stimulation is an important and practical issue. Future work 
will consider these complications, which must also include analysis 
of HF stimulation at frequencies approaching the same time scales 
as cellular dynamics. 

Limitations 

In most clinical settings of interest, HF stimulation is applied in 
the form of an external electrical field stimulus. In this study, we do 
not account for the influence of an external electrical field nor 
account for the spatial extent of the nerve axon. Such levels of 
details are significant for studies of local neural conduction block 
[8-14], as spatial gradients in the extracellular space create virtual 
electrodes resulting in non-uniform HF stimulation through the 
nerve axon [37]. It has been shown in multi-compartment models 
that somatic and axonal firing can become decoupled during 1 00- 
200 Hz HF stimulation, such that somatic quiescence does not 
necessarily preclude activation in neuronal processes [38]. 
Simulation studies in one-dimensional nerve axons have shown 
that as the amplitude of a Idlohertz-range HF stimulation is 
increased, the system can transition from regimes of conduction 
block to rapid firing several times, such that a strict conduction 
block threshold is not clearly defined [39]. As such, non-uniform 
stimulation could lead to conduction block in one region of a 
neuron and rapid activation in another. Further studies of 
kHohertz-range HF stimulation in more spatially-detailed neuronal 
models are necessary to investigate these complex issues. 

The HH model of the giant squid nerve axon is a classical 
model of an excitable cell, highly studied and well-characterized, 
and thus it was a reasonable biopliysically-detailed ionic model to 
characterize the influence of HF stimulation using the method of 
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averaging approach. However, several more detailed neuronal 
models relevant to mammalian physiology have been described, 
incorporating more detailed and multiple sodium, potassium, and 
calcium currents [40-42] . Indeed, the interaction between voltage- 
gated calcium channels and calcium-mediated synaptic transmis- 
sion may be important for understanding the influence of HF 
stimulation on network activity and designing an optimal therapy 
and warrants further study. 

Finally, the simulation results presented here are deterministic 
and do not account for stochastic fluctuations inherent in neuronal 
signaling at both the cellular and subcellular levels. Indeed, studies 
have shown that noise-induced firing can be enhanced by HF 
stimulation for sufficiendy large noise levels, termed vibrational 
resonance [17], closely related to the well-known phenomenon of 
stochastic resonance [43]. Further work is necessary to investigate 
the influence of stochastic fluctuations on spiking in biophysically- 
detailed averaged neuronal models. 

Methods 

Derivation of AFHN model 

We derive the AFHN model equations, following the approach 
in [19]. See [18,19] for a more details. For HF stimuli with large 
frequencies oj » 1 , the period of HF oscillations is much less than 
the characteristic time scales of the FHN neuron. Therefore, we 
seek to eliminate the HF stimulus term poj COS (ojt) from Eq. la 
and obtain an autonomous system which approximates the 
original system on the time scale of the FHN neuron. First, we 
change the variables in Eqs. la and lb, substituting 



K + psin (cot), 



w= W, 



(13a) 



(13b) 



and derive the following equations for V and W: 

V=f[V+psm(o3t)\-W+I, (14a) 



dv _ 1 
dx Inoj 



dw 
dz 



e 

2nw 



{f'[v + p sin (6)]- w + I}de (16a) 



[v + psm(e) + P-yw]de (16b) 



After calculating the integrals and returning to the original time 
scale, the averaged system is as given above in Eqs. 2a and 2b. 
Importantly, we note that the assumption implicit in stating Eq. 
1 3a, specifically that the original system voltage v can be expressed 
as the sum of slow varying V and high frequency term psin(ojt), 
underlies an important conclusion of the method of averaging 
theorem [44]: an equilibrium point in the averaged system (e.g., 
rest or depolarization block) corresponds to a periodic solution in 
the original system (due to the additional HF term superimposed 
on top of the averaged system equilibrium). Similarly, a periodic 
orbit in the averaged system (e.g., repetitive firing) corresponds to 
a more complex oscillation or tori, observed in the v traces of the 
FHN and AFHN models in Figure 1. 

Firing frequency and amplitude in the AFHN model 

Following a similar approach as described in [24], expressions 
for the firing frequency and action potential amplitude in the 
AFHN model are derived as follows. Recall from Eqs. 2a and 2b 
that the v-nuUcline F{v,w) = 0 is cubic in shape. Therefore, over 
the finite range of values for H'e[)Vi,H'2], there are three solutions of 
the equation F(v,M') = 0, which we can denote by v= K_(vP), 
v= Vo(w), and v= K+(m') (see Figure 1 1). The minimal value of vv 
for which K_(vP) exists is H'l, the maximal value of U' for which 
V+(w) exists is Wj, both of which are functions of p and /. 
v=V+(w), the left and right branches, are termed the stable 
branches of the v-nullchne, and v= Vq(w), the middle branch, is 
termed the unstable branch, because, in the limit that v is much 
faster than w, a steady-state located on the (un)stable branch is 
(un)stable. 

The locations of (vi,ivi) and (V2,ii'2) are given by the local 
minimum and maximum, respectively, of the I'-nuUcline, given by 



W = e{V + p sin (cot) + /? - )- W). 



Mathematically, we are interested in the limit gj— » oo, for a fixed 
p. By rescaling time T = cot, we can transform the system to 



dV 



dW 
dx 



--co-^{f[V^p sin {x)\-W + I] 



(15a) 



--(o-'^e{V + psin(x) + fi-yW) (15b) 



and 



WII2-- 



Because v is fast compared to w, v rapidly moves between stable 
branches of the v-nuUcline, V+(w). We can approximate the 
period of an oscillation T by the time required to travel along the 
two stable branches V+{w) 



The variables V and W vary slowly relative to the periodic 
function sinr, due to the small parameter n'^'«l. According to 
the method of averaging [15], an approximate solution to the 
system can be obtained by averaging over the fast periodic 
function, and the averaged variables v and w satisfy the following 
ODEs: 



dt- 



'A 



dt 



(17) 



where points A-D are indiciated in Figure 1 1 . Along the stable 
branches, the dynamics of if are determined by 
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v = V-{w) Vo{w) V+{w) 




P = Ps 1 0 



Figure 11 . Period of a stable limit cycle in the AFHN model. For 

p = 0, the period T of the stable limit cycle is approximately the sum of 
the time required to traverse the two stable branches of the v-nullcline, 
V+{]v), denoted by the points A — B and C — D, respectively. As p 
increases, the amplitude and period of the stable limit cycle, indicated 
by the second set of points labeled A' — D', both decrease. See text for 
description of other variables in figure. Parameters: /= 1.6. 
doi:10.1371/journal.pone.0081402.g011 



^ = I+f'(n){v-n) (22) 

where the dot indicates differentiation with respect to time and 
the prime indicates differentiation with respect to v, such that 

f'(-v,) = c(p)-vl 

and w = Wf Solving Eq. 22 with the initial condition 
v(/ = 0) = v., 

v(t) = -v,-J—{l-txplf'{vM)- (23) 



We set v(?) = ^(p), the v threshold for eliciting an action 
potential, and after rearranging, we arrive at the strength-duration 
curve relationship in Eq. 8. 

To evaluate Eq. 8, we must determine the dependence of the 
threshold C on p. Many studies have discussed the absence of a 
weU-defmed threshold in the FHN model [45-47]. Izhikevich 
notes that canard trajectories following the repelling slow manifold 
provide the best approximation to the excitability threshold [45]. 
Recent work has show that this manifold is well approximated by 
inflection sets (regions of flow lines with zero curvature in the 
phase plane) [46]. For simplicity, we approximate the threshold ^ 
by the middle solution of 



w=G(V+ (w),w) = G+ (w), 
and, therefore, Eq. 17 is equivalendy given by 
^ dw. 



(18) 



(19) 



Note that both terms of the integral are positive, since 
G+(vi')>0 and G_(m')<0 over the range of values [m'i,iv2]. The 
frequency of oscillations is given by \/T. 

From Figure 1 1 , we observe that the amplitude of the limit 
cycle — and thus, of the action potential — is given by — Vc, the 
difference between the values of v at points A and C, respectively. 
Va is the non-repeated root solution of 



f('^A) + I = w\, 



and Vc is similarly the non-repeated root solution of 



(20) 



f(v) + I = w,, 



(24) 



which, as shown in Figure 12 (left panel), does reasonably well- 
approximate the action potential threshold. We are interested in 
the threshold from a non-stimulated state, i.e., 7 = 0. Therefore, 
Eq. 41 can be written as 



1 



— C(p)v + H', 



:(v-V,)(v- + zlv + r), 



(25a) 



(25b) 



where Eq. 25b is implied, since v, by construction is a solution 
of Eq. 25a. Matching coefficients, z) = v, and r= — SvVj/v,, and 
using the quadratic formula, is given by 



-V, — a/v^ + 12vi',/v. 



(26) 



/(Vc) + / = U'2. (21) 



Strength-duration curve in the AFHN model 

An approximation for the strength-duration for the AFHN 
model is derived as follows. Near the steady-state (Vt,Wt), the 
dynamics of v can be approximated by 



where the dependence of ^ on p is due to the dependence of v, 
and vP, on p. Further, from the definition of £,{p), if the terms 
under the radical equal 0, v, is equal to a critical value, which we 
will call — above which the threshold is iU-defmed, that is, the 
system cannot be excited by a brief perturbation from the steady- 
state. Using Eq. 3b, V{ is the real solution of the cubic equation 
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p = 0 




Table 1. Hodgkin-Huxley model current parameters. 



Parameter 


Definition 


Units 


Value 


gNa 


maximum Na^ current conductance 


mS/cm" 


120 


gK 


maximum current conductance 


mS/cm" 


36 


gL 


maximum leal< current conductance 


mS/cm^ 


0.3 


^Na 


Na"'^ current reversal potential 


mV 


115 


Ek 


current reversal potential 


mV 


-12 


El 


leak current reversal potential 


mV 


10.6 


y,a 


resting potential 


mV 


-80 


C 


capacitance 




1 



Parameters for ionic currents in the HH model. 
doi:l 0.1 371 /journal.pone.0081402.t001 



Figure 12. Sub- and super-threshold brief stimuli in the AFHN 
model. (Left) For p = 0, c indicates the threshold for evoking an action 
potential. Two trajectories starting near c are identified by arrows: (1, 
left arrow) when the initial condition v(? = 0)<c, v returns quickly to the 
resting potential v,; (2, right arrow) when v(t = Q)>c, an action 
potential is evoked: the system follows a counterclockwise trajectory, 
quickly approaching the y-nullcline, following the right stable branch 
until reaching the right knee, quickly reaching the left stable branch, 
and then returning to rest. (Right) When p = p,,, C = i'{, the critical value 
above which an action potential cannot be evoked by a brief 
perturbation from the steady-state. 
doi:10.1371/journal.pone.0081402.g012 



yv^ + 12vj + 12/5 = 0. 



(27) 



where v is the shifted transmembrane potential, in which the 
resting potential Vrcsi has been subtracted. The standard 
parameters for the HH model are given in Table 1 . 

Neuronal network model and architecture 

A network of = 1 00 AHH neurons are simulated by adding a 
synaptic current to the AHH model, such that the v dynamics of 
the ;-th neuron are governed by the following equation: 



Note that does not depend on p, only the system parameters. 
Using Eqs. 25a and 27, for all values of p, if 



V/ = [h- gNanr'hivi- EMc)-gKn'(vi- Ek)- gd^^i- El) 



(29) 



and therefore, 



where 



P>Pc = 




the system cannot be excited by a brief stimulus pulse, illustrated 
in Figure 12 (right panel). 

Hodgkin-Huxley model equations 

The equations governing the dynamics of the gating variables 
m,h, and n are given by 

°'"'= exp(2.5-0.1v)-l ^'"='^'^P^-^'/^^^ 
.. = 0.07exp(-v/20)A,= ^^^^^^3-l^ 



and Sex and Si„ are the set of presynaptic neurons with 
connections to neuron i, with excitatory and inhibitory, respec- 
tively, synapses. Sji is the averaged gating variable for the 
postsynaptic conductance, and assumed to be an instantaneous, 
sigmoidal function of the presynaptic cell potential with a 

Table 2. Synaptic current parameters. 



Parameter 


Definition 


Units 


Value 


gsy,i 


maximum synaptic conductance 


mS/cm" 


0.3 


E„ 


excitatory synapse reversal potential 


mV 


80 


Esvn in 


inhibitory synapse reversal potential 


mV 


-12 


V,y„ 


presynaptic cell potential threshold 


mV 


50 


k,„, 


threshold parameter 


mV 


2 



Parameters for excitatory and inhibitory synaptic currents In a network of AHH 
model neurons. 

doi:1 0.1 371 /journal.pone.0081 402.t002 
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threshold V^yn [48], that is 



where 



271 



s.y-_.[vj + psm (e)/C]d6, 



(30) 



where the variance of the time fluctuations of the average 
transmembrane potential, A{t), 



= <(A(0-<A(0»'>, 



^oo(v) = 



1+ exp[-(v- V,y„)/k„.„] 



(31) the variance of the time fluctuations of the individual 

transmembrane potentials v,(?). 



Parameters are in Table 2. The specific synaptic connections 
between neurons are determined randomly, as follows. The 
number of presynaptic connections to the i-th neuron is drawn 
from a Gaussian distribution with mean = 1 0 and standard 
deviation ff = 1 , rounded to the nearest whole number. The 
presynaptic neuron indices j are chosen at random. The type of 
each synapse, excitatory or inhibitory, is determined at random, 
such that the probability of an excitatory synapse is Pexcii£[0,l]- 
Electrical activity is evoked in the neural network by applying a 

200-/iA/cm^, 0.1-ms apphed current in 50 randomly selected 
neurons at time t = 0. 

The collective activity of the neural network can be represented 
by the pseudo-electroencephalogram (pEEG) [49], given by A(;), 



(32) 



the transmembrane potential averaged over all neurons. The 
frequency-domain representation of the pEEG is computed by the 
Fast Fourier Transform. 

The synchrony of the electrical actix-ity in the network is given 
by the synchrony measure ^^[OJl [50], 



N 2^i= 



(33) 



and <x(?)> = (l/r) 



:<(v,(0-<f,(0>r>, 



x(t)dt denotes time-averaging over the 



duration of the simulation T. Note diat if V;(/) are identical for all 
then x=^- 

Numerical simulations 

All numerical simulations were performed in MATLAB. For 
simulations of the AHH model, the modified gating variable rate 
functions (Eqs. 12e-12f) were pre-calculated for a given value of/? 
for vg[- 100,200] mV (K„;G[- 180,120] mV), and values were 
linearly interpolated from look-up tables during simulations. 
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